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ABSTRACT 

The expected gamma-ray flux coming from dark matter annihilation in dwarf spheroidal 
(dSph) galaxies depends on the so-called ‘ J-factor’, the integral of the squared dark matter 
density along the line-of-sight. We examine the degree to which estimates of J are sensitive 
to contamination (by foreground Milky Way stars and stellar streams) of the stellar-kinematic 
samples that are used to infer dark matter densities in ‘ultrafaint’ dSphs. Applying standard 
kinematic analyses to hundreds of mock data sets that include varying levels of contamination, 
we And that mis-classifled contaminants can cause J-factors to be overestimated by orders of 
magnitude. Stellar-kinematic data sets for which we obtain such biased estimates tend 1) to 
include relatively large fractions of stars with ambiguous membership status, and 2 ) to give 
estimates for J that are sensitive to specific choices about how to weight and/or to exclude 
stars with ambiguous status. Comparing publicly-available stellar-kinematic samples for the 
nearby dSphs Reticulum 11 and Segue 1, we And that only the latter displays both of these char¬ 
acteristics. Estimates of Segue I’s J-factor should therefore be regarded with a larger degree 
of caution when planning and interpreting gamma-ray observations. Moreover, robust inter¬ 
pretations regarding dark matter annihilation in dSph galaxies in general will require explicit 
examination of how interlopers might affect the inferred dark matter density profile. 

Key words: astroparticle physics — (cosmology;) dark matter — Galaxy; kinematics and 
dynamics — 7 -rays; general — methods; miscellaneous 


1 INTRODUCTION 

Owing to their proximity, large dynamical mass-to-light ratios and 
low astrophysical backgrounds (Aaronson 1983; Mateo 1998; Mc- 
Connachie 2012), the Milky Way’s dwarf spheroidal (dSph) satel¬ 
lites have become popular targets in ‘indirect’ searches for parti¬ 
cle dark matter (DM) via observations of high-energy photons that 
may be produced in annihilation events (e.g.. Lake 1990; Evans, 
Ferrer & Sarkar 2004; Strigari 2013). In order to translate a given 
photon flux (or non-detection) into constraints (or upper limits) on 
cross-sections for DM annihilation, one must estimate the density 
of DM in the source. This necessity helps to motivate observations 
and dynamical analyses of the tiny stellar populations that, for most 
of these galaxies, represent the only viable tracers of gravitational 
potentials. 

Stellar populations within the Milky Way’s dSph companions 
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are supported against gravity by random motions—i.e., they do not 
rotate. Therefore, dynamical inferences about their gravitational 
potentials must appeal to statistical analyses of distributions of stel¬ 
lar positions and line-of-sight (l.o.s.) velocities. For the most lu¬ 
minous, ‘classical’ dSphs, estimates of DM densities are based on 
kinematic samples containing hundreds to thousands of stars per 
galaxy (Walker, Mateo & Olszewski 2009), giving relatively tight 
constraints (Strigari et al. 2008a; Charbonnier et al. 2011; Cholis 
& Salucci 2012; Geringer-Sameth, Koushiappas & Walker 2015a; 
Bonnivard et al. 2015a). However, many efforts at indirect detection 
focus on the least luminous, ‘ultrafaint’ satellites, for which stellar- 
kinematic sample sizes are necessarily smaller, ranging from single 
to double figures (Martin et al. 2007; Simon & Geha 2007). Despite 
larger uncertainties about their DM halos, the extreme properties 
of ‘ultrafaints’ can give them extraordinary power in constrain¬ 
ing properties of the DM particle (Geringer-Sameth & Koushiap¬ 
pas 2011; Ackermann & Fermi-LAT Collaboration 2011; Geringer- 
Sameth, Koushiappas & Walker 2015b; Ackermann et al. 2014). 

Among the ‘ultrafaint’ dSphs, Segue I (hereafter ‘Seg I’, 
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Belokurov et al. 2007) and the newly-discovered Reticulum II 
(‘Ret IF, Koposov et al. 2015a; Bechtol et al. 2015) have become 
objects of particular interest in the search for gamma-rays from 
DM annihilation. These objects combine the virtues of 1) large in¬ 
ferred DM densities (dynamical masses enclosed within their cen¬ 
tral ~ 30 pc are ~ 10®Mq) and 2) relatively small distances 
{d ~ 30 kpc). As such, these objects offer perhaps the best oppor¬ 
tunities either to detect photons from DM annihilation^ or, in the 
absence of a signal, place stringent upper limits on cross-sections 
for various annihilation channels. 

Given their propagation into constraints on particle physics, it 
is important to examine robustness of dynamical inferences about 
dSph DM densities. In previous work, we have used mock stellar- 
kinematic data sets to examine the sensitivity of dynamical analyses 
to modelling assumptions such as spherical symmetry, parametri- 
sations of DM and stellar density profiles, and behaviour of the 
velocity anisotropy profile (Bonnivard et al. 2015c). In general, we 
find that violations of simplistic modelling assumptions give rise to 
systematic errors that dominate random errors only for the larger 
samples characteristic of ‘classical’ dSphs. For ‘ultrafaints’, such 
systematics are well contained within the relatively large statistical 
errors associated with smaller kinematic samples. 

The conclusions from our previous work hold at least when 
samples are free from contamination by non-member stars, as were 
the mock data sets analysed therein. However, given what is often 
substantial overlap among distributions of observables (e.g., posi¬ 
tion, velocity, chemical abundances), even the largest and most pre¬ 
cise stellar-kinematic data sets do not let us identify with absolute 
certainty which stars belong to the dwarf galaxy and which are con¬ 
taminants in the Galactic foreground. 

Here, we extend our previous work to consider the impact of 
sample contamination on estimates of dSph DM densities and the 
ensuing expectations for indirect detection of DM particles. Specif¬ 
ically, we test the reliability of estimates of the DM annihilation 
factor J — the l.o.s. integral of the squared DM density (Section 2) 
— on a new suite of mock data sets that include varying degrees 
of contamination by populations that mimic Galactic foreground as 
well as tidal streams in which the dSph may be embedded (Sec¬ 
tion 3). Applying analyses that are similar to those used to infer 
DM densities for real dSphs, we examine the severity of systematic 
errors introduced by sample contamination, and we identify char¬ 
acteristic features of data sets that are particularly susceptible to 
such errors (Section 4). Finally, we use these results to gauge the 
reliability of dynamically-inferred J-factors for the specific dSphs 
Seg I and Ret II (Section 5). We stress that this work relies on the 
use of a large number of mock data sets, requiring the analysis to 
be fast. For this reason, the results presented in this study rely on 
a standard Jeans analysis in which the contamination estimation 
and the DM profiles reconstruction are dealt with in two separate 
steps. An alternative Bayesian analysis, using a global likelihood 
which allows to constrain simultaneously the two, is presented in 
Appendix A. However, the relative computational expense of this 
method inhibits the extensive testing to which we have subjected 
the more standard analyses employed in this and previous work. 


^ Based on public data from the Fermi Large Area Telescope, Geringer- 
Sameth et al. (2015) report an excess of emission in the direction of Ret II 
at energies ~ 2 — 10 GeV with significance > 2.3(7. Based on proprietary 
data, the Drlica-Wagner et al. (2015) finds an excess at the same energies 
in Ret II, but assigns the signal significance of just 0.3 (t. Revisiting the 
public Fermi-LAT data. Hooper & Linden (2015) confirm the detection of 
Geringer-Sameth et al. (2015) and assign it similar significance. 


In the few cases where we have tested it and can compare di¬ 
rectly against the performance of the standard analyses, we find no 
significant improvement (e.g., in heavily-contaminated mock data 
sets where the standard analyses give strongly biased results for J- 
factors, we find similarly strong biases using the global likelihood). 
As summarized in Section 6 , this indicates that a two step analysis 
is sufficient for assessing contamination impact on J-factors. Im¬ 
proving the algorithms used to separate members from contaminant 
in dSphs, and/or getting more data should be the priority to obtain 
more robust constraints on contaminated objects such as Segue I. 


2 ESTIMATION OF DARK MATTER DENSITY AND THE 
J-FACTOR 

The differential 7 -ray flux received on Earth in solid angle Afl 
from DM annihilation is related to the DM particle mass, the 
thermally-averaged annihilation cross section {(Jv), the energy 
spectrum of annihilation products, and the quantity 

J(AH) = j I Pd]v[(/, H) d/dlJ . (1) 

J An J 

The ‘J-factor’, calculated here with the public code CLUMPY^ 
(Charbonnier, Combet & Maurin 2012; Bonnivard et al. 2016), is 
the integral of the squared DM density along l.o.s. I and over solid 
angle AH = 27r x [1 — cos aint], with aint the angle over which 
the 7 -ray signal is observed. For a given particle physics model and 
measurement of 7 -ray flux, the ensuing constraints on DM proper¬ 
ties depend directly on J. For dSph galaxies, the J-factor is best 
constrained at the critical integration angle Qc = 2 x r^/d, where 
rh is the half-light radius and d is the distance to the dSph (Walker 
et al. 2011; Charbonnier et al. 2011). 


2.1 Jeans analysis 


Following many previous analyses (e.g., Strigari et al. 2007; Essig 
et al. 2010; Charbonnier et al. 2011), here we estimate DM density 
profiles, and hence J-factors, by modelling stellar-kinematic data 
under assumptions of spherical symmetry, dynamical equilibrium 
and negligible rotational support. The Jeans equation then relates 
observables (stellar positions and velocities) to the gravitational po¬ 
tential: (Binney & Tremaine 2008) 


1 d 
V dr 


(uv^) + 2 


/3ani(r)tl^ 

r 


GM{r) 


( 2 ) 


Here, v{r) is the stellar number density profile, v^ir) is the radial 
velocity dispersion profile, /3ani(p) = 1 — v'g/v'^ is the anisotropy 
profile of the stellar velocity dispersion, and M(r) the mass— 
assumed to be dominated by DM—enclosed within radius r. Cur¬ 
rent observations resolve neither internal proper motions nor rela¬ 
tive distances of stars in a given dSph, confining empirical infor¬ 
mation to three dimensions of phase space (projected position and 
l.o.s. velocity). After solving Eq. (2) and projecting along the l.o.s., 
the (squared) stellar velocity dispersion at projected radius R reads 
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with I(R) the projected stellar number density profile. 

Given a stellar-kinematic data set, we compare the l.o.s. stellar 


^ http://lpsc.in2p3.fr/clumpy 
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velocities to a model velocity dispersion, ap{R), computed using 
parametric forms for the velocity anisotropy /5ani (f) and DM den¬ 
sity profile pDM(f) (Bonnivard et al. 2015c). For the DM density 
profile, pDM(f), we adopt the Einasto profile: 

PDM(f) = p- 2 exp|-- ( —) -1 

where a controls the logarithmic slope, r _2 is a scale radius and 
p _2 is a scale density. 

We consider real data for known dSphs as well as mock data 
for which the true DM density profile is known. Before fitting the 
kinematic data for real dSphs, we estimate I(R) by fitting a model 
of the form 

= ,5, 

JR Vr —R 

—i.e., the sum of the l.o.s. projection of i/(r) and a uniform back¬ 
ground density—to the unbinned distribution of positions of stars 
identified as red giant candidates in photometric survey data. Fol¬ 
lowing Bonnivard et al. (2015c, see their Section 2 for details), we 
model the stellar density profile as a broken power law: 



where a* controls the sharpness of transition from inner log-slope 
7 * to outer log-slope /3* at scale radius rj. When analysing mock 
data sets, we adopt the true profiles for v{r), 1{R), and /?ani(f), 
in order to isolate systematic errors arising from contamination of 
stellar-kinematic samples. For the two real ‘ultrafaint’ dSphs con¬ 
sidered in Section 5, we assume constant anisotropy, /3ani = Po, but 
have confirmed that the results would be similar if we used more 
flexible models for /3ani (e.g. Baes & van Flese profile (Baes & van 
Flese 2007), see Bonnivard et al. 2015a). 

We estimate the free parameters using either a Markov Chain 
Monte Carlo (MCMC) engine, or a minimisation algorithm with 
bootstrap resampling (see Section 2.3). Thus, our fits have three 
free parameters (p_ 2 , f _2 and a) when applied to mock data, and 
four (now including a constant anisotropy /3o) when applied to real 
data. For all parameters we adopt the same priors as Bonnivard et al. 
(2015a, see their Table 1). 


2.2 MCMC analysis: Jmedian and credible intervals 

Following Bonnivard et al. (2015a,b), we use the GreAT MCMC 
engine (Putze & Derome 2014) to efficiently sample the multi¬ 
dimensional parameter space. This procedure returns posterior 
probability density functions (PDFs) for model parameters. Via Eq. 
(1), we can compute median values and credible intervals (CIs) for 
J-factors at any integration angle (e.g., Charbonnier et al. 2011). 

For our MCMC analysis, we adopt a likelihood function that 
includes stellar velocity measurements as discrete data points (Stri- 
gari et al. 2008b; Martinez et al. 2009)®. Rather than discarding 
observations from stars suspected of being contaminants, we can 
include all data and weight the observation by a probability 
of membership. Pi, that can be estimated from observables under 
specific assumptions about member and contaminant populations 

^ An alternative method would be to bin the velocites in order to estimate 
the observed velocity dispersion profile, crahs{R)^ and compare it to the 
model a-p{R) (Walker et al. 2009a; Charbonnier et al. 2011). Using large 
sets of mock data, Bonnivard et al. (2015a) found that the two methods give 
similar results. 


(e.g., Klimentowski et al. 2007; Wojtak & Lokas 2007; Walker et al. 
2009b; Martinez et al. 2009). We adopt a likelihood function that 
has natural logarithm 


log£ = 


JV„ta,rs .. JVstars 

- ^ P,log(2^)-- ^ P.log(a®(7?.) + A®J 
■ 1=1 2 = 1 

1 Pi{Vi — (b))® 
"2 ^ a^{R,) + AlS 


i.e. we assume that the distribution of l.o.s. stellar velocities v is 
Gaussian, centred on the mean stellar velocity (v), and that the dis¬ 
persion of velocities at radius R comes from both the intrinsic dis¬ 
persion o-p{R) and the measurement uncertainty A„. 


2.3 Bootstrap analysis: ( J)bootstrap and (^bootstrap 

We also use in this work a bootstrap estimator (Efron 1982) for 
the J-factor and its statistical uncertainty. For a given kinematic 
data sample, we generate N — 500 resamples by drawing, with 
replacement, n stars from the n observations of the original data 
set. For each resample, we find the best-fitting model using a x® 
minimisation algorithm, with the function 

x" = -21og£, (8) 

log/1 being the log-likelihood from Eq. (A5). We then compute 
the J-factor of each resample. The mean and dispersion over the 
N = 500 values are the bootstrap estimators of the J-factor 
((J)bootstrap) and its statistical uncertainty (abootstrap). 

The motivation for this approach is that it can reveal sensitivity 
to the presence of a few outliers in cases where sample contamina¬ 
tion is important. The estimated J-factor might change drastically 
from one realization to the next, depending on whether or not the 
contaminants are drawn, and this can result in larger statistical un¬ 
certainties (than the associated MCMC-derived uncertainties) on 
the J-factor. 


3 MOCK DATA DESCRIPTION 

In order to investigate the impact on estimates of J-factors due to 
contamination of stellar-kinematic samples by non-member stars, 
we generate thousands of mock data sets that each sample mixtures 
of three simulated stellar populations tracing a gravitational poten¬ 
tial dominated by DM. 


3.1 Simulated populations 

The first simulated stellar population represents bona fide mem¬ 
bers of an ‘ultrafaint’ (resp. ‘classical’)-like dwarf galaxy^. For 
this population, following the procedure described by Charbonnier 
et al. (2011), we draw phase-space coordinates randomly from dis¬ 
tribution functions that satisfy the collisionless Boltzmann equa¬ 
tion (Binney & Tremaine 2008, see Section 3.2.1 for details). The 
second population represents contamination from a tidal stream in 
which the ‘ultrafaint’ may be embedded, as may be the case if the 
object formed as the satellite of a more massive and less dense 
‘classical’ dwarf that was more easily disrupted by Galactic tides 


The ‘classical’-like mock dSphs are studied in Bonnivard et al. (2015a) 
(see their Appendix C), and will not be further discussed in this paper. 
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(Belokurov et al. 2009). The third population represents contami¬ 
nant stars in the Galactic foreground, which often have the same 
colours and magnitudes as the red giants targeted in spectroscopic 
surveys of dSph galaxies. 

In order to generate a given mock stellar-kinematic data set, 
we first determine the contribution from each stellar population hy 
drawing three random numbers. The first, drawn from the interval 
logiQ Nmin ^ logio ^ ^ logio -^max with Uniform probability, 
sets the overall sample size, N. For realisations corresponding to 
‘ultrafaint’ (resp. ‘classical’) dSphs, we adopt log A^min = 1-5 and 
logA^max = 2.0 (resp. logA^min = 2.5, logA^max = 3.5). The 
second random number, drawn from the interval 0.1 ^ /i ^ 0.9 
with uniform probability, gives the fraction of N stars drawn from 
either the member or stream population. The third, drawn from the 
interval 0.1 ^ f 2 ^ 0.9 with uniform probability, gives the fraction 
of those fiN stars that are drawn from the member population. 
Thus the numbers of stars drawn from member, stream and Galactic 
foreground populations are A'^/1/2, Nfi{l — /2) and A^(l — /i), 
respectively. 

3.2 Observables 

Then, for each simulated star corresponding to a given population, 
we draw observables (position, velocity, metallicity, surface grav¬ 
ity) from plausible distributions for that population. 

3.2.1 Positions and velocities for members 

As mentioned above, for dwarf members we randomly draw pro¬ 
jected positions and l.o.s. velocities from the same phase-space dis¬ 
tribution functions that we have previously used to generate mock 
data sets in previous studies (Charbonnier et al. 2011; Bonnivard 
et al. 2015c). We have chosen four models in order to capture 
a range of behaviours corresponding to various plausible proper¬ 
ties of dSph DM halos and the phase-space distributions of their 
stellar tracers. One has a large (~ 100 pc) ‘core’ of constant 
central DM density; two have centrally ‘cusped’ density profiles 
(d log p/d log r = — 1 as r —> 0 ), as in halos produced in cosmo¬ 
logical N-body simulations (Springel et al. 2008); the fourth has 
a shallower cusp (dlogp/dlogr = — 0.6 as r —>■ 0 ), mimicking 
a possible effect of stellar feedback in low-mass halos (Brooks & 
Zolotov 2014). One model has a stellar velocity dispersion profile 
that increases with radius; another’s velocity dispersion profile de¬ 
creases with radius, and two have velocity dispersion profiles that 
are approximately constant. Among the four models, global stellar 
velocity dispersions range from 4.5 km s “ ^ to 11 km s “ ^. We draw 
1000 mock data sets for each model, giving a total of 4000 mock 
‘ultrafaints’. 


3.2.2 Positions and velocities for contaminants 

In order to allow for varying degrees of overlap with the Galactic 
velocity distribution (held constant; see third bullet point below), 
we shift the mean l.o.s. velocity of dwarf members by a value that 
we draw randomly from the interval —200 ^ (v) ^ 4-200 km s“^. 

• For stream and Galactic foreground stars, we draw projected 
positions from distributions that are uniform over the field sub¬ 
tended by dwarf members. 

• For stream stars, we draw l.o.s. velocities from a Gaussian dis¬ 
tribution that has the same mean velocity as the dwarf members. 


but—in order to represent a more massive parent system—twice 
the velocity dispersion (Belokurov et al. 2009). 

• For Galactic foreground stars, we draw l.o.s. velocities from 
the distribution that is expected to be observed along the l.o.s. to 
the Fornax dSph, based on the Galactic model of Besancon (Robin 
et al. 2003). 

All mock observables are scattered by 4 km s“ ^ in velocity in order 
to mimic realistic observational errors (e.g., Simon & Geha 2007). 

3.2.3 Metallicity and surface gravity 

For the member stars within a given realisation, we adopt a mean 
metallicity that is drawn from a distribution that is uniform between 
—3.5 ^ ([Fe/H]) ^ —1.5. For a given member star within that 
realisation, we then draw metallicity from a Gaussian distribution 
with the adopted mean and dispersion of cr[pe/H] = 0.5 dex, typical 
of observed ultrafaints (Kirby et al. 2010). 

• In all realisations, we draw metallicities for stream stars from a 
Gaussian distribution with mean {[Fe/H]) = —1.5 and dispersion 
'^[Fe/H] = 1 dex. 

• For member and stream stars, we draw surface gravity from 
a Gaussian distribution with mean (logj^Q g) — 1.5 (cgs units) and 
dispersion criog 9 = 1 dex. 

• For Galactic foreground stars, we draw metallicities from 
a Gaussian distribution with mean —0.5 dex and dispersion 1.5 
dex, and surface gravities from a Gaussian distribution with mean 
(logio 5 ) = 4 and dispersion 1 dex. 

All mock observables are scattered by 0.2 dex in metallicity and 
0.2 dex in surface gravity in order to mimic realistic observational 
errors (e.g., Kirby et al. 2010). 

3.2.4 Membership probabilities 

For each star in a given mock data set, we estimate probability of 
membership using the expectation-maximisation (EM) algorithm 
described by Walker et al. (2009b). The EM algorithm iteratively 
uses estimates of parameters (e.g., means and dispersions) that 
specify distributions of observables (e.g., velocity, metallicity, sur¬ 
face gravity) in order to evaluate membership probabilities for in¬ 
dividual stars, then uses those probabilities to update estimates of 
model parameters. 

Because we use a single Gaussian distribution to represent 
the velocity distribution of member stars, our implementation of 
the EM algorithm effectively assumes that the velocity disper¬ 
sion profile is flat, introducing potential biases when the dSph has 
a radially-variable velocity dispersion profile. Despite our use of 
four different dynamical models whose velocity dispersion profiles 
range from decreasing to flat to increasing (Section 3.2.1), we did 
not find any systematic errors that were specific to one particu¬ 
lar kind of model. We conclude, therefore, that our current study 
specifically examines the effect of sample contamination and does 
not depend strongly on the nature of the underlying dynamical 
model. 


4 J-FACTORS FOR CONTAMINATED MOCK DATA 

4.1 Uncertain membership status and contamination 

We first look at the impact of contamination on the estimation of 
the membership probabilities with the EM algorithm. Based on 
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Figure 1. Top: fraction /inter of plausible members (membership proba¬ 
bility P > 10“®) whose membership status is ambiguous (0.05 < P < 
0.95) for our 4000 mock data sets, as a function of the fractions of the sam¬ 
ple contributed by contamination from the Milky Way {x axis) and stream 
(y axis). Bottom: accuracy of estimated J-factors (MCMC/cut-95% analy¬ 
sis, aint = ac), for the 545 mock dSphs with /inter > 0.1. The colour 
scale shows deviation of the estimated J-factor from the true value. The 
sizes of the symbols correspond to the size of the 68% Cl on the estimated 
J-factor. 


their study of 27 mock data sets of different levels of contamina¬ 
tion from simulated foreground populations, Walker et al. (2009b) 
conclude that the algorithm reliably separates contaminants from 
dSph members, except when samples are small (Wstars < 30) and 
suffer heavy contamination (i.e., when there is a small number of 
members buried among a large sample of contaminants). 

Here, using our 4000 ‘ultrafaint’-like mock dSphs, we focus 
on the impact of stars that show ambiguous membership status, i.e., 
have intermediate membership probabilities. We define a new di¬ 
agnostic, /inter, as the fraction of plausible members that have am¬ 
biguous status. Specifically, /inter is the ratio of the number of stars 
with 0.05 < P < 0.95 to the number of stars with P > 10“®. 
This quantity is useful because it can be computed for real data 
sets. For our 4000 mock data sets. Figure 1 shows how /inter de¬ 
pends on intrinsic levels of contamination by simulated Milky Way 
and stream stars, where the fractions of stream and Milky Way 
stars in a given mock data set are given by /mw = 1 ~ /i and 
/stream = /i(l “ /a), respectively (see Section 3.1). 

Strong levels of contamination, particularly from the Milky 
Way, tend to result in large fractions of stars having ambiguous 
membership status. Therefore, the quantity /inter can serve as a 
proxy for contamination levels in real dSphs, and will be discussed 
for the cases of the ‘ultrafaint’ dSphs Ret II and Seg I in Section 5. 
Henceforth, in order to focus on situations where sample contam¬ 
ination is important, we consider only the 545 mock data sets for 
which /inter > 10% (see bottom panel of Fig. 1). 


4.2 Impact of contamination on J-factors 

Four different procedures For each of the 545 mock data sets 
with /inter > 10%, we perform both the MCMC and bootstrap 
analyses described in Sections 2.2 and 2.3. For both analyses we 
handle the estimated membership probabilities in two different 
ways, both of which have been employed by various authors in 
analyses of real dSph data. In the first case, we use all data points 
and weight each by its probability of membership according to Eq. 
(A5); henceforth, we refer to this procedure as ‘P-weighted’. In 
the second, we discard all observations for which the probability 
of membership fails to exceed a fiducial threshold of P > Pthresh, 
and we give equal weight to all observations that exceed this thresh¬ 
old. In the second case, we still use Eq. (A5) as the log-likelihood, 
but only after reassigning membership probabilities to P = 0 
for all stars originally below the threshold and to P = 1 for all 
stars originally above it. Following several published analyses (e.g., 
Charbonnier et al. 2011), we adopt Pthresh = 0.95, and hence¬ 
forth refer to this procedure as ‘cut-95%’. Note that the member¬ 
ship probabilities themselves are likely to carry large uncertainties, 
but the latter are not estimated by the different reconstruction algo¬ 
rithms (e.g., the EM algorithm). We are therefore left to deal with 
their point estimates. 

Considering the MCMC and bootstrap analyses and both sce¬ 
narios for handling membership, we obtain four independent esti¬ 
mations of the J-factor for each mock ‘ultrafaint’, which we com¬ 
pute here at aint = Oc = 2 x r^/d (Walker et al. 2011; Charbon¬ 
nier et al. 2011), with the half-light radius and d the distance to 
the mock dSph®. We then compare the estimated J-factors to the 
‘true’ values that we calculate directly from the DM density profiles 
that we used to generate the mock data. In what follows, all statis¬ 
tical uncertainties are la, representing 68% CIs for the MCMC 
analyses and (Tbootstrap for the bootstrap analyses. 

J-factor overestimation The bottom panel of Figure 1 depicts ac¬ 
curacy of the J-factor estimates in the /mw — /stream plane for the 
545 mock dSphs. These particular results correspond to the cut- 
95% procedure with MCMC parameter estimation (see below for 
a comparison of results from all four procedures). Dark blue sym¬ 
bols correspond to cases where J is underestimated; these are rare 
and tend to have relatively large CIs. Colours ranging from green 
to red represent cases where J is overestimated by factors as large 
as ~ 10®. These cases tend to occur when samples have a large 
degree of contamination, near the line /mw + /stream = 1. The 
largest overestimates are driven by foreground contaminant stars 
that are mistaken for members, with the EM algorithm returning 
P > 0.95. The velocity dispersion of the Milky Way stars is sig¬ 
nificantly larger than the one from the member stars, and the mis- 
classification of even a few such contaminants can mimic large 
DM densities and yield large J-factors. Alarmingly, even when J 
is severely overestimated, the CIs can be relatively small. Since 
real dSphs with large and tightly-constrained J-factors receive the 
greatest weight in subsequent inferences of cross-sections for DM 
annihilation, it is important to recognise when such estimates may 
be biased by sample contamination. 

The top panel of Eigure 2 shows distributions of 
logic= ctc for the 545 mock data sets 
with /inter > 10%, and for the four different analysis procedures. 
All four yield similar distributions, with long tails extending toward 

® Note that our conclusions do not depend on the distance to the mock 
dSph. 
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Figure 2. Top\ distributions of logjQ( J/Jtrue) at for the 545 

mock dSphs, each analysed with four different procedures that differ ac¬ 
cording to method of parameter estimation (MCMC vs. bootstrap) and treat¬ 
ment of membership probabilities (weighting of individual data points vs. 
discarding low-probability members). All procedures yield long tails cor¬ 
responding to overestimates of J. Bottom', distributions of relative differ¬ 
ences between the four procedures; the ‘reference’ procedure is taken to be 
the case in which we estimate parameters via MCMC and discard stars for 
which estimated probabilities of membership are smaller than 0.95. 


overestimated J, reaching values as large as logjg(J/Jtrue) ~ 5. 
However, the most extreme cases of overestimation are due to 
strong levels of Milky Way contamination, which give very large 
velocity dispersions that would be easily identified when analysing 
real dSphs. If we select for instance the models for which the 
mean stellar velocity dispersion Op is smaller than 20 km s“^, the 
maximum overestimation reduces to ~ 3 orders of magnitude. The 
distribution remains however skewed toward overestimated values 
because of the contamination from both Milky Way and stream 
stars, which is not the case with contamination-free mock data 
(Bonnivard et al. 2015c). 

Comparing results from different procedures In order to exam¬ 
ine sensitivity to choice of procedure, we systematically compare 
results from the four different procedures that we applied to each 
mock data set. For ease of comparison, we refer the result from 
each procedure to the one obtained from our MCMC analysis of the 
cut-95 % procedure that was used to construct Figure 1. The bottom 
panel of Figure 2 show distributions of logj^g) J/obtained for 
our 545 samples. 

For most of the models, the four procedures give similar 
results, such that distributions of logjQ(J/peak near zero. 
However, when procedures give different estimates of J, we tend 
to get larger values when using the membership probabilities as 
weights. This is true for both MCMC and bootstrap methods (red 


and orange distributions respectively). The reason is that lower- 
probability members tend to populate tails of the members’ veloc¬ 
ity distribution. Analyses that give non-zero weight to these stars 
will necessarily operate on broader velocity distributions and infer 
larger DM masses. 

Figure 3 provides a concrete example, showing results from 
our analysis of one particular mock data set. The bottom panel 
shows how velocity and membership probability depend on pro¬ 
jected distance from the dSph center. The fraction of plausible 
members with ambiguous status is /inter = 0.35. Depending on 
how we treat these stars, the velocity dispersion profile changes dra¬ 
matically: it is smaller by a factor ~ 2 when removing the stars with 
P < 0.95 (middle panel).® The top panel shows how our estimate 
of J depends strongly on how we treat these stars. The P-weighted 
procedure gives non-zero weight to the ambiguous stars, which in 
this case are non-members. As a result, the P-weighted procedure 
overestimates J by two orders of magnitude and, despite this inac¬ 
curacy, has relatively small CIs (dashed red curves and filled orange 
zone, for MCMC and bootstrap analyses respectively). In contrast, 
the cut-95% procedure gives zero weight to ambiguous cases and 
recovers the true value of J within its broader CIs. 

Lessons from analysis of mock data Our experimentation with 
mock data sets teaches us that when stellar-kinematic sam¬ 
ples are small and heavily contaminated, J-factors tend to be 
overestimated—sometimes by orders of magnitude—as the result 
of mis-classification of even a few non-member stars. This conclu¬ 
sion stands, even if a more sophisticated method is used to simul¬ 
taneously estimate the membership status and constrain the DM 
profile. This method is presented in App. A. Being more computa- 
tionnaly demanding, it could not be tested extensively on the mock 
data, and such detailed analysis is beyond the scope of this paper. 
It was however found to give similar error bars despite propagat¬ 
ing correctly the membership uncertainties to the DM profile esti¬ 
mation — the uncertainties remain dominated by the sample size 
— and was not able to better retrieve the correct J-factor. For the 
pathological H0234 mock model discussed above, it performs even 
worse than the cruder standard analysis, overshooting the correct 
J-factor by several orders of magnitude. For these reasons, the two- 
step approach (extracting membership probabilities and then using 
them in the Jeans analysis) seems adequate for identifying cases 
where contamination is likely to be particularly problematic. 

In the context of this approach, cases of severe overestimation 
tend to exhibit characteristics that can be identified in real data sets. 
Warning signs include: 

• relatively large fractions of stars with ambiguous membership 
status, as quantified by /inter; 

• sensitivity of J-factor estimates to the exact procedure chosen 
for addressing membership (e.g., using membership probabilities 
for all stars vs. discarding observations below some threshold). 


5 SEGUE I AND RETICULUM II 

Having gained intuition from our systematic analysis of mock data 
sets, we now consider real stellar-kinematic data sets that are pub¬ 
licly available for two ‘ultrafaint’ dSphs: Seg I (Simon et al. 2011) 

® We also show the best-fit models obtained in each case from the MCMC 
analysis. As we employ an unbinned analysis, these fits are shown only for 
illustration purpose. 
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Figure 3. Top: J-factor as a function of integration angle for the mock 
ultrafaint dSph H0234. The true J-factor, in solid black, is compared to es¬ 
timates obtained with four different analyses that differ based on choice of 
procedure for parameter estimation (MCMC vs. bootstrap) and treatment 
of membership probabilities (weighting of all observations by membership 
probability vs. discarding all observations with P ^ 0.95). Filled regions 
and dashed curves correspond to la uncertainties obtained via bootstrap 
and MCMC procedures, respectively. Middle: line-of-sight velocity dis¬ 
persion profile for H0234, estimated either using membership probabili¬ 
ties as weights (red) or considering only observations for which probabil¬ 
ity of membership is > 0.95 (black). The best-fit models from the un¬ 
binned MCMC analyses are shown for illustration purpose. Bottom: for 
the same mock dSph, deviation from the mean velocity, |u — (v)\ (blue 
to red colours), as a function of membership probability, P (only stars with 
P > 10“® are plotted), and projected radius, R. 


and Ret II (Walker et al. 2015; Simon et al. 2015; Koposov et al. 
2015b). 


5.1 Reticulum II 

We first focus on the kinematic sample of Ret II. Following Bon- 
nivard et al. (2015b), we use the kinematic data of Walker et al. 
(2015), which contains estimates of velocities, effective tempera¬ 


tures, surface gravities and metallicities for 37 red giant candidates 
along the l.o.s. to Ret II'. We quantify membership probabilities 
using the EM algorithm of Walker et al. (2009b). The sum of mem¬ 
bership probabilities is ~ 16.7, where 16 stars have P > 0.99 and 
another has P « 0.7. 

The bottom panel of Figure 4 shows deviations from the mean 
velocity, \v— (v) \ (blue to red colours), as a function of membership 
probability P and projected radius R, for stars with P > 10“^. 
Only the star with P « 0.7 has ambiguous membership probabil¬ 
ity (due to its relatively large metallicity), and no obvious velocity 
outlier is visible. The velocity dispersion profiles obtained either 
by weighting the contribution of each star by its value of P, or by 
selecting only the stars with P ^ 0.95, are very similar (middle 
panel of Figure 4). Applying our four different procedures gives 
confidence that our estimate of J for Ret II is robust against con¬ 
tamination, as the P-weighted and cut-95 % analyses show good 
agreement for both MCMC and bootstrap analyses (top panel of 
Figure 4). 


5.2 Segue I 

For Seg I, we analyse the stellar kinematic data set of Simon et al. 
(2011, ‘Sir hereafter), who measure velocity and equivalent width 
of the calcium-triplet absorption feature for 393 stars in the direc¬ 
tion of Seg I. For the ~ 180 of these stars that have colours and 
magnitudes consistent with Seg I membership, Sll estimate two 
different membership probabilities, using either the EM algorithm 
or a Bayesian technique described by Martinez et al. (2011). For 
straightforward comparison with our analyses of mock data sets, 
we first focus on the membership probabilities obtained from the 
EM algorithm. Following Sll, we remove from the sample two 
photometric variable stars (one of which Sll identify as a velocity 
variable), three additional stars that Sll identify as velocity vari¬ 
ables, and one additional velocity outlier that Sll identify as having 
a strong effect on estimates of velocity dispersion. 

Even after removing these stars, the remaining sample for 
Seg I shares characteristics with the mock data sets for which we 
found estimates of J to be unreliable. The fraction of stars with am¬ 
biguous membership probability is /inter = 0.19 (bottom panel of 
Figure 5). Moreover, for both MCMC and bootstrap procedures for 
parameter estimation, our estimates of J for Seg I depend strongly 
on whether we use the cut-95 % or the P-weighted procedure to 
handle membership probabilities (top panel of Figure 5). For ex¬ 
ample, using the P-weighted procedure, we obtain a J-factor at 
Qint = 1° that is more than two orders of magnitude larger than 
the result from the cut-95% procedure. This behaviour is strikingly 
similar to that of the mock data set illustrated in Figure 3, and is due 
to the increase of velocity dispersion at large radius when including 
the ambiguous stars (middle panel of Figure 5). 


5.3 A deeper look at Segue I 

Our different procedures return significantly different estimates of 
J for Seg I, indicating a lack of robustness against contamination 
of the stellar-kinematic sample. Let us examine this sensitivity in 
greater detail. 


^ We confirm that we obtain the same results if we use instead the stellar- 
kinematic samples of Simon et al. (2015) or Koposov et al. (2015b). 
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Figure 4. Same as Figure 3, except for the real ‘ultrafaint’ dwaif Ret II, 
using the stellar-kinematic data of Walker et al. 2015). Only one star has 
ambiguous membership probability (P Ri 0.7), and all four procedures for 
estimating J give consistent results, indicating robustness against contami¬ 
nation of the stellar-kinematic sample. 

5.3.1 Bayesian algorithm for membership probabilities 

First, we repeat our analysis of SH’s data for Seg I, but using the 
membership probabilities that Sll obtain via the Bayesian analy¬ 
sis described by Martinez et al. (2011). This gives a larger number 
of stars with ambiguous membership probability (/inter increases 
from 0.19 to 0.46; bottom panel of Figure 6). The velocity disper¬ 
sion profiles obtained with the different procedures are very dif¬ 
ferent (middle panel of Figure 6): the velocity dispersion at large 
radius, obtained with the cut P ^ 0.95, is three times smaller than 
when using the P values as weights. Accordingly, the differences 
between J-factors from our four procedures become even larger 
than what we obtained previously (top panel of Figure 6): the P- 
weighted analysis gives a value of J that is more than three orders 
of magnitude larger than the one obtained using the cut-95 % anal¬ 
ysis. Interestingly, the cut-95% bootstrap analysis yields extremely 
large uncertainties (green area in the top panel of Figure 6), indi¬ 
cating that the estimation of J may be unduly influenced by the 
presence of a few outliers in the sample. 


Figure 5. Same as Figure 4, but for Seg I, using tbe stellar-kinematic sample 
of Simon et al. (2011). 


5.3.2 Sensitivity to membership probability threshold 

Our choice of threshold Pthresh = 0.95 for selecting high- 
probability members is arbitrary. We now examine the impact 
of our choice of Pthresh on the estimation of Seg Fs J-factor, 
still using our MCMC procedure and Sll’s Bayesian estimates of 
membership probabilities. Figure 7 shows how results change as 
we vary the threshold from Pthresh = 0.1 to 0.95. Thresholds 
of Pthresh Ss 0.5 give Small J-factors with large uncertainties. 
Thresolds of Pthresh < 0.5 give larger J-factors with much smaller 
uncertainties, the natural consequence of including more stars with 
large deviations from the mean velocity. 

Sll themselves note that their estimate of Seg Fs velocity 
dispersion is sensitive to the inclusion/exclusion of one particu¬ 
lar star, SDSSJ100704.35-1-160459.4, which is a 6cr velocity out¬ 
lier but, since it is close to the center of Seg I, receives P > 0.95 
from the EM algorithm and P ~ 0.5 from the Bayesian algorithm. 
Following their lead, we excluded this star from all of our analy¬ 
ses (see Section 5.2), so it cannot be responsible for the sensitiv¬ 
ity we find with regard to choice of procedure. Instead, we find 
strong sensitivity to whether we include or exclude another star. 
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Figure 6. Same as in Figure 5, but with the Bayesian estimation of the 
membership probabilities (Simon et al. 2011; Martinez et al. 2011). The 
black square in the bottom panel highlights a star whose exclusion changes 
drastically the J-factor reconstruction. See text for discussion. 


SDSSJ100659.95+160408.7® (yellow dot highlighted by a black 
square in the bottom panel of Figure 6), for which the Bayesian 
algorithm gives P ~ 0.39. For cases with Pthresh < 0.4, Fig¬ 
ure 7 illustrates this sensitivity by showing the effect of removing 
this one star from the sample that is otherwise specified by Pthresh- 
Open black circles show results from analyses that include the star, 
and yellow circles with rings show results from analyses that ex¬ 
clude it. The difference can be as large as ~ 2 orders of magnitude 
in J. 


5.3.3 Understanding discrepancies among previously derived 
J-factors 

Several published analyses of Seg I have used S11’s data to estimate 
J-factors that are larger and less uncertain than we have recently 
reported (Bonnivard et al. 2015a, ‘B15’ hereafter). In that work, we 
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18 
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Figure 7. J(aint = cuc = 0.14°) and 95% CIs from P-weighted (red 
cross) and P-unweighted MCMC analyses as a function of Pthresh’ the 
threshold for which stars with P < Pthresh ^re discarded from the anal¬ 
ysis. The difference between filled yellow circles and black circles is the 
removal of a single star shown with a similar symbol in the bottom panel of 
Figure 6 (star at P = 5.5 pc and with P = 0.39). 
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used exclusively the ‘cut-95%’ procedure with MCMC estimation 
of parameters. In order to understand apparent discrepancies with 
respect to BIS’s result for Seg I, logJ(0.25°) = I 6 . 8 I 2 j [GeV^ 
cm“®], we have reproduced two published analyses of Seg I which 
use similar techniques. 

First, Essig et al. (2010, ‘ElO’ hereafter) apply a similar Jeans 
analysis to the 66 stars from S11 for which the EM algorithm gives 
P > 0.8, all of which receiving an equal weight. ElO assume an 
Einasto profile for the DM density, an isotropic velocity dispersion 
(Pani ~ 0) and a Plummer light profile ((a, /3, 7 ) = (2,5, 0) in Eq. 
6 ). Compared to our setup, their prior ranges on the DM density are 
narrower, with e.g. a G [0.14 : 0.3] while we use a € [0.12 : 1] 
(see Bonnivard et al. 2015a). ElOreport log J(0.25°) = 19.1±0.6, 
with CIs based on a fit of a Gaussian function to the posterior 
PDF that they obtain for J. Relaxing our membership threshold to 
Pthresh = 0.8 and adopting the same assumptions as ElO, we can 
reproduce their estimate so long as we also follow their procedure. 
However, we find that the PDF for J is significantly non-Gaussian, 
with long tails that extend toward small values. While a Gaussian 
fit gives log J(0.25°) = 18.8 ± 0.5, consistent with ElO’s result, 
the median value (with CIs enclosing the central 68 % of probabil¬ 
ity) of the actual PDF is log J(0.25°) = 18.leading to much 
larger uncertainties. Thus the discrepancy between B 15’s and ElO’s 
estimates of Seg I’s J-factor is due to 1) their use of a less restric¬ 
tive member sample, 2) their fit of a Gaussian function to a non- 
Gaussian PDF, and 3) their more restrained ranges on DM density 
priors. 

We have also repeated the analysis of Geringer-Sameth, 
Koushiappas & Walker (2015a, ‘G-S15’ hereafter), who use the 
S11 sample and enforce a membership threshold Pthresh >0.5 (G- 
S15’s paper erroneously reports using a threshold of 0.95 for Seg I; 
this stricter threshold was used only for their analysis of ‘classical’ 
dSphs). We can reproduce G-S15’s estimate of J if we adopt their 
modelling assumptions and use Pthresh = 0.5. We find that the dis¬ 
crepancy between G-S15’s estimate of log J(0.5°) = lO.SG^'lg'gg 
and B15’s is due to 1) different thresholds for membership prob¬ 
abilities, and 2 ) their choice of anisotropy prior, which favours 
/3ani = 1 (radial anisotropy), while we employ a flat prior® on /lani 
on the range [—9 : 1]. 


® This star is not among the seven cases of ambiguous membership that 
S11 discuss in their Appendix. 


® Bonnivard et al. (2015c) found that the exact parametrisation of /Jani has 
no strong impact on the reconstruction of the J-factor of ‘ultrafaint’ dSphs, 
but the shape of the prior can be important. 
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Figure 8. Comparison of the J-factors of Seg I, obtained with the several 
tests of this paper, to the values for the closest dSphs galaxies from Bon- 
nivard et al. (2015a,b). The J-factor can vary from ~ 10^® to ~ 10^^ 
GeV^ cm“® from one analysis to another. On the other hand, Ret IFs J- 
factor was found to be robust against the different tests. 


6 SUMMARY 

In order to help guide efforts at indirect detection of DM parti¬ 
cles in dSph galaxies, we have examined sensitivity of expected 
gamma-ray fluxes to contamination of the stellar-kinematic sam¬ 
ples on which estimates of DM densities are based. Using a large 
suite of mock stellar-kinematic data sets that include varying de¬ 
grees of contamination by Galactic foreground and stream popula¬ 
tions, we compare estimated to true values of the ‘J-factor’, which 
is proportional to the expected gamma-ray fluxes from DM anni¬ 
hilation. We find that, for the small stellar-kinematic samples that 
are available for ‘ultrafaint’ dSph galaxies, contamination by fore¬ 
ground stars can cause J-factors to be overestimated by several 
orders of magnitude. 

In order to gauge the level at which contamination affects es¬ 
timates of J for real dSphs, we have identified characteristics of 
mock data sets for which estimates of the J-factor are unreliable. 
Such data sets tend 1) to have a large fraction of stars with ambigu¬ 
ous membership status, as quantified by estimates of membership 
probabilities, and 2) to give different results depending on details 
of arbitrary choices about whether to include or exclude specific 
stars with ambiguous membership status. 

It is difficult to judge what fraction of the known dSph galax¬ 
ies are vulnerable to these sorts of complications. For the 8 known 
‘classical’ dwarfs we studied in Bonnivard et al. (2015a), we find 
that only Fornax shows specific signs of contamination. We have 
considered the stellar-kinematic data sets that are available for two 
of the ‘ultrafaint’ dSph galaxies that have been the focus of several 
indirect detection investigations: Ret II and Seg I. The available 
samples for Ret II suggest a relatively small degree of contamina¬ 
tion and contain relatively few stars with ambiguous membership 
status; as a consequence, our estimates of J for Ret II do not exhibit 
sensitivity to various possible procedures for addressing member¬ 
ship. In contrast, the available sample for Seg I has a larger frac¬ 
tion of foreground contaminants, and a substantially larger fraction 
of stars with ambiguous membership status. Furthermore, we find 
that estimates of J for this object can vary by several orders of 
magnitude depending on the inclusion/exclusion of a few stars with 
ambiguous membership status. In general, future analyses of dSph 
galaxies should check explicitly for such sensitivity to interlopers. 

Figure 8 summarises the sensitivity to contamination that we 
find when analysing Seg I’s kinematic sample in various ways, and 


compares the range of estimated J-factors to those we obtain for 
other nearby dSphs (Bonnivard et al. 2015a,b). For procedures that 
give non-zero weight to marginal members in the wings of Seg I’s 
velocity distribution, we obtain J-factors that are larger and have 
relatively small CIs compared to other dSphs; similar results from 
similar analyses have motivated deep gamma-ray observations of 
Seg I (e.g., Aliu et al. 2012; Aleksic et al. 2014). On the other hand, 
procedures that use stricter membership criteria give estimates for 
J that are significantly smaller and more uncertain than obtained 
for other dSph galaxies^®. Depending on choice of procedure, for 
Seg I we can recover estimates of J-factors spanning ~ 3 orders of 
magnitude, covering the range of previously published values. The 
results of the more sophisticated, full Bayesian analysis (App. A) 
gives a high J value, but we note that this method appears to suffer 
from the same (and sometimes worse) drawbacks than the simpler 
analysis in presence of strong contamination levels. We emphasize 
however that a more detailed testing of this approach, beyond the 
scope of this paper, would be necessary to understand the origin 
of these drawbacks. We conclude that estimates of J-factors for 
Seg I should be regarded with extreme caution when planning and 
interpreting indirect detection experiments. 

Further spectroscopic measurements of Seg I could help re¬ 
solve this issue, by improving the quality and the amount of data. 
This task, as well as improving algorithms for separating members 
from any kind of foregrounds, should receive a high priority given 
the importance of this object. 
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APPENDIX A: ALTERNATIVE PROCEDURE EOR 
ASSESSING MEMBERSHIP AND ESTIMATING J 
SIMULTANEOUSLY 

Both the ’cut-95’ and ’P-weighted’ analyses suffer from the facts 
that 1) stellar membership is evaluated under a different model (i.e., 
one that assumes the velocity dispersion of members is indepen¬ 
dent of position) than the one that is used subsequently to infer 
p(r), and 2) uncertainties in membership probabilities are not taken 
into account. Therefore we also consider an alternative analysis that 
allows for simultaneous evaluation of membership and inferences 
about p(r) (and hence J(aint))- Specifically, we consider a mixture 
model under which P{R, V, Z\9) represents the joint probability 
density of stellar position R, line-of-sight velocity V and metallic- 
ity Z, given a model specified by the vector 6 of free parameters. If 
a star is drawn randomly from a mixture of Np populations, where 


For Ursa Major II (UMa II), Willman I (Wil I) and Coma, we could not 
repeat our tests for sensitivity to choice of procedure for handling member¬ 
ship probabilities, as the relevant stellar-kinematic samples are not publicly 
available. 
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(unobserved) variable M € {1, 2,..., Np } indicates that star’s pop¬ 
ulation, then 

Np 

P{R, V, Z\e) = ^ P{R, K z, M\e) 

M=l 

Np 

= ^ p{M\e)P{R,v,z\M,e) 

M=1 

Np 

= ^ P(M|e)P(_R|M,0)P(V,Z|M,_R,e). (Al) 

M=1 

In the first line, we marginalise over the unknown variable M. In 
each of the last two lines we apply the chain rule of probability: 
P{A, B) = P(A|P)P(P). Here, P{M\6) is the probability that 
the star is drawn from population M, given the model. P{R\M, 0) 
is the probability density of position, given the model and that the 
star is drawn from population M. P{V, Z\M, R,0) is the joint 
probability density of velocity and metallicity, given the model and 
that the star is drawn from population M and has position R. 

The probability density P{R\M, 9) can be written in terms of 
projected stellar density profile, Im{R), for the population spec¬ 
ified by M. Assuming spherical symmetry, a circle of radius R 


encloses Nm{R) = 2n 
as seen in projection, and 


P(RIM,9) = — 


S Im {S) dS of this population’s stars 


d 

\Nm{R)] 

d 

pTT ffsiM(S) dNl 

dR 

_Nm{oo) 

dR 

AAf(oo) 


2tvRIm{R) 

Nm{oo) 


(A2) 


The probability P{M\9) is the fraction, /ai, of stars that belong to 
population M: 


-/Vtot 


(A3) 


where Atot = Ai(oo)-|-A2(oo)-|-...-|-Ap(oo) is the total number 
of stars, including all populations. Substituting Eqs. (A2) and (A3) 
into Eq. (Al), we obtain 




n ( E (RdPiVj ,Z,\R„M, 9) \. (A5) 


P{R, V, Z\9) = ^ V Im{R)P{V, Z\R, M, 9). (A4) 

M=1 

Then a data set consisting of independent observations of po¬ 
sition, velocity and metallicity for N stars has likelihood 

JL/'n -TTR. 

jC = 

-*■ -*■ \ /Vt„* 

M=1 

In order to analyse stellar-kinematic data sets for dwarf spheroidals 
in the presence of foreground contamination, here we assume 
Ap = 2, allowing for member (M = 1) and foreground (M = 2) 
populations. We model the 3D stellar density using Eq. (6). We 
assume the stellar density of foreground stars is uniform over the 
relatively small field subtended by the member population, such 
that the foreground population has projected density I 2 , a constant. 
Then the fraction of observed stars that are members is 

4^fcir?r(^)r(5^) 


4^ P 2 

/i = — rvi{r)dr = 


aiN 


(A6) 


where r(a:) is the gamma function. The fraction of observed stars 
that are foreground is 


/2 = l-/l = --, 


(A7) 


where i?max is the maximum radius of the (circular) field. 

Finally, we assume that the velocities and metallicities of both 
member and foreground populations follow bivariate normal distri¬ 
butions. 


PiV,Z\M,R,9)=Af2{pM,Z;M), 


with mean vectors 


Um = 


({V)m\ 


and covariance matrices 



(A) + 5^ 
0 


0 




(A8) 


(A9) 


(AlO) 


whose diagonal elements are broadened by observational errors 5v 
and &z (we implicitly assume that velocity and metallicity are un¬ 
correlated). We assume the velocity dispersion of the foreground 
population is uniform across the field of view, but we allow the ve¬ 
locity dispersion of the member population to vary with position. 
Specifically, we use Eq. (4) to model p(r), and then use the Jeans 
equation to compute the variance of the projected velocity distribu¬ 
tion as described in Section 2.2. 

Our mixture model is fully specified by 18 free parameters: 


/logio[pJ(Af0/pc3)]\ 

logio[fo/pc] 

a 

d 

7 

-logroll - 
fi 

logic [fi/pc] 
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(l")l 

{V}2 

{Z)i 

{Z)2 

<^V,2 

O’!,! 

V <^ 1.2 


(All) 
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